Dúvidas pra ver com o Guilherme

Dúvidas em azul eu acho que já entendi, mas deixei aqui a resposta para que a Aisha do futuro possa consultar também. Dúvidas em vermelho ainda precisam de resposta.

Dúvidas de tradução:

Resposta: * u.m.c. milhão (é a unidade que o Bacen usa para algumas variáveis) * variação % mensal (tipo quando pega o PIB e calcula o de hoje menos o de ontem dividido pelo de ontem) * IPCA - price index ou inflation rate ? * BM&FBOVESPA S.A. - Bolsa de Valores, Mercadorias e Futuros (BM&FBOVESPA) * Capital income * Labor income * Capital-Labor ratio * Interest rate (o que fazer com SELIC, deixo SELIC mesmo?) * 3 month treasury bill rate * Price index ou Inflation rate? * Real Effective Interest Rate Index (Índice da Taxa de Câmbio Efetiva Real) * Precisa traduzir as fontes dos dados? Por exemplo, Banco Central departamento de estatísticas ou a Bovespa?

Dúvida 2: Queria discutir sobre raizes unitárias…

Como o modelo é TVP, por construção ele não “acredita” em uma distribuição estacionária?

Resposta: Isso também. Mas o fato de usar um passeio aleatório como lei de movimento dos coeficientes naturalmente impõe uma estrutura não estacionária no modelo.

1. Dúvida 2: Tá, mas aí como que a gente faz com a função impulso resposta? Pelo que eu entendo, a função impulso resposta usa os coeficientes calculados (no caso, estou usando a média) em um determinado $t$ para dar o choque em uma variável e ver o que ocorre com as demais, mas ele não muda os coeficientes… Como que fica isso? Eu fui procurar na documentação do pacote que estou usando como que é calculada a FIR. Tem três formas de fazer (eu estou usando a terceira, que é a mesma maneira do Primiceri e Del Negro (2015)). O modelo do Primideri para uma série temporal n-dimensional é dado por:

\[y_t = c_t + B_{1,t}\cdot y_{t-1} + \ldots + B_{p,t} \cdot y_{t-p} + A_t^{-1} \cdot \Sigma_t\cdot \varepsilon_t\] (a especificação completa pode ser vista aqui).

No primeiro cenário, a FIR é calculada considerando o vetor $u_t \equiv A_t^{-1} \cdot \Sigma_t\cdot \varepsilon_t$. A função estima qual o efeito de um choque de um desvio padrão em um dos elementos $u_t$ (i.e. $u_{t.j}=1$ e $u_{t,i}=0$, para todo $i \neq j$) afeta $y_{t+h}$, em comparação com $u_t$ sendo um vetor de zeros. Esta abordagem não é considerada realista pois usualmente $A_t^{-1} \cdot \Sigma_t$ não será uma matriz diagonal (os choques são correlacionados). Nesta abordagem, tanto os elementos de $c_t$ como de $B_{j,t}$ são mantidos fixos nos seus valores estimados para $t$;
No segundo cenário, a FIR estima o impacto de uma unidade em algum dos elementos de $\varepsilon_t$. Os elementos de $c_t$, $B_{j,t}$, $A_t$ e $\Sigma_t$ são mantidos fixos nos seus valores estimados para $t$. Esta abordagem depende da ordem das variáveis.
O cenário 3 é igual ao cenário dois mas os elementos de $\Sigma_t$ são substituídos pela média dos valores em todos os períodos. Esta forma é a que foi utilizada por Primiceri em 2005 e Primiceri e Del Negro em 2015.

Na prática o que o algoritmo faz é calcular a FIR para cada uma das estimativas obtidas pelo algoritmo MCMC e então calcula a média e os quantis desses valores.

Resposta:

Dúvida 3: Sobre o X-13 ARIMA SEATS. Tem uma opção de usar uma transformação nos dados antes de ajustar o modelo SARIMA. O que o procedimento faz é uma transformação Box-Cox. Mas daí a minha dúvida é o seguinte: em geral o modelo testa se deve fazer o log versus fazer nada (o default é esse). Só que a série sem sazonalidade, quando plotada contra a série original está sempre na mesma escala (e no caso do PIB, que tem tendência, o bicho continua com tendência). Essa transformação de log é o quê, então?

Eu li e não achei resposta pra isso:

Cap 3 do DeJong,
http://documentation.sas.com/?docsetId=etsug&docsetTarget=etsug_x13_syntax20.htm&docsetVersion=14.2&locale=en,
https://www.imf.org/external/pubs/ft/qna/pdf/2017/chapter7.pdf (ver nota de rodapé 14 na página 10)
https://www.census.gov/ts/papers/gettingstartedx13acc_winx13.pdf (ver página 6)
ftp://ftp.ibge.gov.br/Contas_Nacionais/Contas_Nacionais_Trimestrais/Ajuste_Sazonal/X13_NasContasTrimestrais.pdf
https://usefulr.wordpress.com/2014/02/17/adjusting-chinese-new-year-effects-in-r-is-easy/

O pouco que eu acho que entendi é que a função por default usa um critério (baseado na verossimilhança) para definir se o modelo SARIMA ajustado é multiplicativo ou aditivo. Mas isso ainda não me ajuda a entender porque a escala da variável dessazonalizada não é em log, então alguma coisa ainda está confusa…

Resposta:

Dúvida 4: Quando mudei a taxa de câmbio para taxa real efetiva, as funções impulso resposta do Câmbio sobre a razão capital trabalho passaram a ser significativas. Porém, está com sentido de que “aumentos na taxa de câmbio aumentam a razão capital trabalho”. Só que é o mesmo formato que a taxa de juros tem, “choque positivo nos juros, aumenta a razão capital-trabalho”. O problema está em eu esperar que juros e câmbio tenham relação inversa, isto é, eu esperaria que um aumento na taxa de juros fizesse o câmbio diminuir - se o juros aumenta fica mais interessante para investimentos estrangeiros e com isso há um aumento de moeda estrangeira fazendo com que a taxa de câmbio fique menor - então se juros alto $\Rightarrow$ câmbio baixa, eu esperaria que a relação que cada um deles tem com a razão capital trabalho fosse invertida.

Por outro lado, se a taxa de câmbio aumenta, os investimentos em moeda estrangeiro ficam mais rentáveis e os rendimentos do capital ficariam maiores, pensando exclusivamente em capital financeiro. Mas se pensar em termos de capital físico, se uma pessoa precisa importar uma máquina de sorvete (de ANDRADE, C.R., 2018) e o câmbio aumentou, ela vai precisar gastar mais dinheiro para comprar a máquina. Por outro lado, quem exporta, vai ganhar mais. Argh. ~~Merda de economia aberta~~.

Resposta:

Dúvida 5: Eu estou com problema para interpretar o coeficiente do PIB defasado na equação dele, pois está negativo. Isso quer dizer que quanto maior a taxa de variação em $t-1$, menor será a taxa de variação em $t$?

Resposta:

Dúvida 6: Na documentação do pacote, na parte de previsão, ele calcula algo que ele chama de média condicional aos parâmetros ($\mu_{T+j}^{(\ell)}$). Isso seria como a média incondicional do modelo naquele instante? Tô meio confusa, porque o modelo não é estacionário… [C:/Users/aisha/Pictures/hmm.jpg]

Resposta:

Dúvida 7: Com relação ao MAPE e MSE, eu calculo para cada variável individualmente, correto? Não é tipo uma medida “global”…?

Resposta:

Introduction

This code has all the calculations used in the work “The Impact of Monetary Policy on Income Inequality: Evidences from Brazilian Economy”. It is divided in three parts: the first one deals with data preparation and some descriptive statistics; second part estimates a VAR model and third part assess the convergence diagnostics for the MCMC algorithm.

To estimate the VAR model, I used the bvarsv do Fabian Kruger package. This two files are very helpful: bvarsv: An R implementation of the Primiceri (2005) model for macroeconomic time series and Replication of figures in Del Negro and Primiceri (2015). For the diagnostics, I used the coda package.

Downloading packages

If mirror 10 (“UFRJ”) is offline, try changing ind = 10 to ind = 1. Do the same thing if you get an error installing the BETS package.

chooseCRANmirror(graphics = FALSE, ind = 10)
if (!require("pacman")) install.packages("pacman")
pacman::p_load(devtools, ggplot2, forecast, BETS, seasonal, seasonalview, bvarsv, lubridate, zoo, stargazer, gridExtra, reshape2, ggfortify, RColorBrewer, scales, quantmod, PerformanceAnalytics, strucchange, coda, foreach, doParallel, knitr, grid, ggpubr, gdata, vars, urca, compiler, DescTools, kableExtra, knitr)

These are some auxiliary functions used later in the code (I took them from bvarsv documentation):

## Some graphical parameters
matplot2 <- function(...){
  matplot(..., type = 'l', lty = 1, lwd = 2, bty = "n", ylab = "")
}

# Stat.helper calculates que quantiles and organizes them in quantile 1, mean, quantile 3 given a vector z
stat.helper <- function(z) c(mean(z), quantile(z, c(0.16, 0.84)))[c(2,1,3)]

# Color palette
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

cols1 <- cbPalette[c(2,4,2)]
cols2 <- cbPalette[c(2,4,6)]

# This one is for IRF analysis
stat.helper2 <- function(z) c(mean(z), quantile(z, c(0.05, 0.25, 0.75, 0.95)))[c(2, 3, 1, 4, 5)]

Data handling

Downloading data

All data was obtained from Brazil’s Central Bank repository (Sistema Gerenciador de Séries do Banco Central - SGS). The following series were used (obs: beginning and end of the series corresponds to the available data, not actual data used in the work and some data were manipulated, for example interest rate were converted to annual - mannipulations are described in code and in the paper appendix):

Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do trabalho (Labor Income)
- Beginning of the series: January, 1992.
- End of the series: February, 2018.
- Unity: u.m.c. milhões.
- Series number in SGS: 7620.
- Source: Brazil’s Central Bank - Statistics Department (BCB-DSTAT).
Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do capital (Capital Income)
- Beginning of the series: January, 1992.
- End of the series: February, 2018.
- Unity: u.m.c. milhões.
- Series number in SGS: 7621.
- Source: Brazil’s Central Bank - Statistics Department (BCB-DSTAT).
Taxa de Juros - Selic (Interest Rate)
- Beginning of the series: July, 1986.
- End of the series: April, 2018.
- Unity: % monthly.
- Series number in SGS: 4390.
- Source: Brazil’s Central Bank - Open Market Operations Department (BCB-Demab).
Taxa referencial de swaps DI pré-fixada (BM&F) - Prazo de 90 dias (fim de período) (3-month treasury bill rate)
- Beginning of the series: September, 1999.
- End of the series: March, 2018.
- Unity: % annual.
- Series number in SGS: 7818.
- Source: BM&FBOVESPA S.A. - Bolsa de Valores, Mercadorias e Futuros (BM&FBOVESPA).
Índice nacional de preços ao consumidor-amplo (IPCA) (Price index)
- Beginning of the series: January, 1980.
- End of the series: March, 2018.
- Unity: Percent monthly variation.
- Series number in SGS: 433
- Source: Brazilian Institute of Geography and Statistics (IBGE).
Monthly GDP
- Beginning of the series: January, 1990.
- End of the series: February, 2018.
- Unity: Brazilian Real (R$) - millions.
- Series number in SGS: 4380
- Source: Brazil’s Central Bank - Economic Department (BCB-Depec).
Interest Rate - USD Dollar (sell) - End of Period - Monthly - ANTIGO
- Beginning of the series: January, 1953.
- End of the series: March, 2018.
- Unity: u.m.c./US$.
- Series number in SGS: 3696.
- Source: Brazil’s Central Bank - Sisbacen PTAX800.
Real Effective Interest Rate Index - Monthly
- Beginning of the series: January, 1988.
- End of the series: March, 2018.
- Unity: Index (100 = Jun/1994).
- Series number in SGS: 11752
- Source: Brazil’s Central Bank - Statistics Department (BCB-DSTAT).

Obs: Para ver o ADF, usei isso aqui: https://stats.stackexchange.com/questions/24072/interpreting-rs-ur-df-dickey-fuller-unit-root-test-results

Downloading the data

# Auxiliary variables, so I don't need to bother when something changes
inicio <- "1996-01-01"
fim    <- "2018-02-28"
inicio_deflator <- "1993-12-31"
inicio_swap     <- "1999-09-01"

# Don't mess with this code
inicio_cambio <- paste(seq(as.Date(inicio), length = 2, by = "-1 month")[2]) # 1 month before the beginning of the other series
inicio_ipca2   <- paste(seq(as.Date(inicio), length = 12, by = "-1 month")[12]) # 12 months before the beginning of the other series
fim_ipca      <- paste(seq(as.Date(fim), length = 2, by = "1 month")[2]) # 1 month after the end of the other series
inicio_swap2 <- paste(seq(as.Date(inicio_swap), length = 2, by = "-1 month")[2]) # 1 1 month before the beginning of the other series

## Obs: In case you have a SQL problem, rollback this two packages and mannually download BETS from CRAN and reinstall it as well
## See: https://stackoverflow.com/questions/43073782/rmysql-system-error-10060 and https://github.com/pedrocostaferreira/BETS/issues/29
# require(devtools)
# install_version("DBI", version = "0.5", repos = "http://cran.us.r-project.org")
# install_version("RMySQL", version = "0.10.9", repos = "http://cran.us.r-project.org")

################################
##### Capital labor ratio ######
################################

trabalho <- BETS.get("7620", from = inicio, to = fim)
capital  <- BETS.get("7621", from = inicio, to = fim)
capital_trabalho <- capital/trabalho

########################
##### Price Index ######
########################

ipca_raw <- BETS.get("433", from = inicio_deflator, to = fim_ipca) 
ipca2 <- ipca_raw
ipca2[1] <- 100

# The following calculation was made following IPEA methodology: http://www.ipeadata.gov.br/iframe_transformacao.aspx?width=1474&height=701
# Transforms in index
for (i in 2:length(ipca2)){
  ipca2[i] <- round(ipca2[(i-1)]*(1+ipca2[(i)]/100),2)
}

# Geometric Average
for (i in 1:(length(ipca2)-1)){
  ipca2[i] <- sqrt(ipca2[i]*ipca2[(i+1)])
}

# Change base to mar/2018
for (i in 1:length(ipca2)) {
  ipca2[i] <- tail(ipca2, n=2)[1]/ipca2[i]
}

# Elimitates data outside the period we are studying
ipca2 <- ts(ipca2,  start = c(1993, 12, 31), frequency = 12)
ano_ipca <- as.numeric(substr(inicio_cambio, start = 1, stop = 4))
mes_ipca <- as.numeric(substr(inicio_cambio, start = 6, stop = 7))
ipca2 <- window(ipca2, c(ano_ipca, mes_ipca))

ano_ipca <- as.numeric(substr(inicio_ipca2, start = 1, stop = 4))
mes_ipca <- as.numeric(substr(inicio_ipca2, start = 6, stop = 7))
ipca_raw <- window(ipca_raw, c(ano_ipca, mes_ipca))

# Removes last observation
ipca_raw <- ipca_raw[1:(length(ipca_raw)-1)]

ipca_acum <- ipca_raw/100 + 1
ipca <- vector()

final <- length(ipca_acum)
for (i in 12:final){
  ipca[(i-11)] <- (prod(ipca_acum[(i-11):i])-1)*100
}

ano <- as.numeric(substr(inicio, start = 1, stop = 4))
mes <- as.numeric(substr(inicio, start = 6, stop = 7))
dia <- as.numeric(substr(inicio, start = 9, stop = 10))

ipca <- ts(ipca,  start = c(ano, mes, dia), frequency = 12) # Date format YYYY MM DD

################
##### GDP ######
################

PIBmensal <- BETS.get("4380", from = inicio_cambio, to = fim)
PIBmensal <- PIBmensal * ipca2 

##########################
##### Exchange rate ######
##########################

# Antigo câmbio nominal
#cambio_raw <- BETS.get("3696", from = inicio_cambio, to = fim)
cambio_raw <- BETS.get("11752", from = inicio, to = fim) # Since this is already an index, I don't need to take the log-diff

###########################
##### Interest rates ######
###########################

selic_4390 <- BETS.get("4390", from = inicio, to = fim)
# Transforms into anual rate
selic <- ((1+selic_4390/100)^(12)-1)*100

swap90 <- BETS.get("7818", from = inicio_swap, to = fim)

Visualizing data

Treating for seasonality and structural breaks

The output of the seasonality tests will appear. The analysis is hidden inside the code.

1. Capital-Labor Ratio

# Calculates the structural breaks
bp_ts <- breakpoints(capital_trabalho ~ 1)
ci_ts <- confint(bp_ts)

# Monthplot aggregates days and months and takes averages
meses <- ggmonthplot(capital_trabalho) + 
  theme_bw() +
  #scale_y_continuous(name = "Capital-Labor Ratio (% - average)") +
  scale_y_continuous(name = "Razão capital-trabalho (%)") +
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))
  #labs(title = "Média diária e mensal da razão capital-trabalho", x = "Mês", subtitle = "Nenhuma transformação")
  #labs(title = "Daily and monthly averages of Capital-Labor ratio", x = "Month", subtitle = "No data transformation")

# Seasonplot has the averages for month, each line represents a year
season <- ggseasonplot(capital_trabalho, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Razão capital-trabalho (%)") +
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))
  #scale_y_continuous(name = "Capital-Labor Ratio (% - average)") +
  #labs(title = "Monthly averages of Capital-Labor ratio", x = "Months", subtitle = "No data transformation")

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Season_CT.pdf", width = 6, height = 3)
season

#dev.off()

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Month_CT.pdf", width = 6, height = 3)
meses

#dev.off()


Data <- seq(as.Date("1996-01-01"), length = length(capital_trabalho), by = "1 month")
df1 <- melt(data.frame(Data,capital_trabalho), id.vars = "Data")
names(df1) <- c("Data", "Variavel", "Valor")

# Graph with structural break (dashed) and CI (red)
p1 <- ggplot(df1[which(df1$Variavel == "capital_trabalho"),], aes(Data, Valor, colour = Variavel)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name = "Razão capital-trabalho (%)") +
        #scale_y_continuous(name="Capital-Labor ratio") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+ 
        geom_vline(xintercept = Data[bp_ts$breakpoints], linetype="longdash")+
        geom_segment(aes(x = Data[ci_ts$confint[1]], y = min(capital_trabalho), xend = Data[ci_ts$confint[3]], yend = min(capital_trabalho)))+
        theme_bw()
p1 <- p1 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), legend.position = "none")

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Quebra_CT.pdf", width = 6, height = 3)
p1

#dev.off()

# We divide the series in pieces (accordingly to the structural break)
capital_trabalho_a <- capital_trabalho[1:bp_ts$breakpoints]
capital_trabalho_b <- capital_trabalho[(bp_ts$breakpoints+1):length(capital_trabalho)]
capital_trabalho_a <- ts(capital_trabalho_a,  start = c(1996, 01, 01), frequency = 12)
capital_trabalho_b <- ts(capital_trabalho_b,  start = c(2004, 02, 01), frequency = 12)

m <- seas(x = capital_trabalho_a, transform.function = "none", regression.aictest = NULL)
#m <- seas(x = capital_trabalho_a, regression.aictest = NULL)

## First diagnostics is using qs() function to check for seasonality after the adjustment
qs(m)

##                    qs   p-val
## qsori        34.12733 0.00000
## qsorievadj   91.53477 0.00000
## qsrsd         0.00000 1.00000
## qssadj        0.00000 1.00000
## qssadjevadj   0.00000 1.00000
## qsirr         0.00078 0.99961
## qsirrevadj    0.00000 1.00000
## qssori       33.13775 0.00000
## qssorievadj  85.86391 0.00000
## qssrsd        0.00000 1.00000
## qsssadj       0.00000 1.00000
## qsssadjevadj  0.00000 1.00000
## qssirr        0.05069 0.97497
## qssirrevadj   0.00000 1.00000

# Click here to see the how the qs statistics is calculated: https://stats.stackexchange.com/questions/148573/the-results-and-specifics-from-the-qs-function-in-r
# The first 7 results are related to the whole series, the last 7 comprehends only the last 8 years
# qsori is the original series
# qsorievadj is the original series corrected for extreme values
# qsrsd is the residual's series of the ARIMA model
# qssadj is the series with seasonal adjustment 
# qssadjevadj is the series with seasonal adjustment corrected for extreme values
# qsirr is the series of the irregular component
# qsirrevadj is the series of the irregular component corrected for extreme values
# by the output, we can conclude that there is no evidence of seasonal components in the adjusted series, the residuals and the irregular component, even when we look just the last 8 years.

# Next step is to diagnose the pre-adjust and the ARIMA model
summary(m)

## 
## Call:
## seas(x = capital_trabalho_a, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AO1996.Jan        -0.35561    0.10377  -3.427 0.000611 ***
## AO1998.Jan         0.67046    0.08634   7.766 8.12e-15 ***
## LS1998.Jul         0.28775    0.05339   5.390 7.06e-08 ***
## AO1998.Aug         0.29547    0.08564   3.450 0.000561 ***
## AO1999.Feb         0.68086    0.09704   7.016 2.28e-12 ***
## AO1999.Mar         0.37531    0.09644   3.892 9.96e-05 ***
## AO2002.Oct         0.49217    0.08564   5.747 9.10e-09 ***
## AR-Nonseasonal-01  0.64477    0.07100   9.081  < 2e-16 ***
## MA-Seasonal-12     0.99965    0.08975  11.138  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (1 0 0)(0 1 1)  Obs.: 97  Transform: none
## AICc: -109.8, BIC: -88.36  QS (no seasonality in final):    0  
## Box-Ljung (no autocorr.): 20.02   Shapiro (normality): 0.9765 .

# We have an level-shift outlier in July, 1998 (Asian crisis?)
# We also have additive outliers in Jan96, Jan98, Aug98, Feb99, Mar99, Oct02
# The estimated model is an ARIMA (1 0 0)(0 1 1) with the non-seasonal AR significative and the MA seasonal significative.
capital_trabalho_2a <- final(m)

m <- seas(x = capital_trabalho_b, transform.function = "none", regression.aictest = NULL)

#m <- seas(x = capital_trabalho_b, regression.aictest = NULL) # se tomar o log, fica melhor
qs(m)

##                     qs   p-val
## qsori        242.32197 0.00000
## qsorievadj   287.80760 0.00000
## qsrsd          0.00000 1.00000
## qssadj        13.25076 0.00133
## qssadjevadj    0.00000 1.00000
## qsirr          6.90372 0.03169
## qsirrevadj     0.00000 1.00000
## qssori       134.89609 0.00000
## qssorievadj  155.37794 0.00000
## qssrsd         4.35854 0.11312
## qsssadj        5.14745 0.07625
## qsssadjevadj   0.00000 1.00000
## qssirr         7.59824 0.02239
## qssirrevadj    0.00000 1.00000

summary(m)

## 
## Call:
## seas(x = capital_trabalho_b, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AO2004.Jun        -0.74155    0.05278 -14.049  < 2e-16 ***
## LS2004.Nov        -0.25272    0.03722  -6.789 1.13e-11 ***
## AO2004.Dec        -0.46336    0.04768  -9.718  < 2e-16 ***
## AO2005.Jan         0.21001    0.05129   4.095 4.22e-05 ***
## AO2005.Jun         0.39668    0.05002   7.930 2.19e-15 ***
## AO2006.Jan         0.41973    0.04679   8.969  < 2e-16 ***
## AO2006.Jun         0.52133    0.04725  11.034  < 2e-16 ***
## AO2006.Dec        -0.19085    0.04160  -4.588 4.47e-06 ***
## AO2007.Jan         0.35202    0.04432   7.943 1.98e-15 ***
## AO2007.Jun         0.36303    0.04432   8.191 2.59e-16 ***
## AO2010.Jun        -0.20396    0.04037  -5.053 4.36e-07 ***
## AO2011.Oct         0.22042    0.03978   5.541 3.01e-08 ***
## AO2011.Dec         0.19854    0.03984   4.983 6.26e-07 ***
## AO2013.Jun        -0.40527    0.04064  -9.972  < 2e-16 ***
## AO2014.Jun        -0.31716    0.04068  -7.797 6.35e-15 ***
## AO2016.Dec         0.24956    0.04116   6.063 1.34e-09 ***
## MA-Nonseasonal-01  0.68849    0.05725  12.027  < 2e-16 ***
## MA-Seasonal-12     0.63001    0.06814   9.246  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (0 1 1)(0 1 1)  Obs.: 169  Transform: none
## AICc: -456.5, BIC: -404.1  QS (no seasonality in final):13.25 **
## Box-Ljung (no autocorr.): 26.23   Shapiro (normality): 0.9833 *

capital_trabalho_2b <- final(m)

###############
# Juntando novamente as duas séries
###############

comb <- ts.union(capital_trabalho_2a, capital_trabalho_2b)
capital_trabalho_final <- pmin(comb[,1], comb[,2], na.rm = TRUE)

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho_final), by = "1 month"), capital_trabalho, capital_trabalho_final)
names(df) <- c("Data", "Original", "Filtrada")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = capital_trabalho, colour = "Original"), linetype = "dashed") +
  geom_line(aes(y = capital_trabalho_final, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Ano", labels = date_format("%Y"))+
  scale_y_continuous(name = "Razão Capital/Trabalho (%)")+
  labs(color="Série") +
  scale_colour_manual("Variável", breaks = c("Original", "Filtrada"), values = c(cores[1], cores[2]))+
  theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7), legend.position = "top")

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-AntesDepois_CT.pdf", width = 6, height = 3)
p1

#dev.off()



# Essa parte ficou de fora
# Now I'm going to test for unit roots and outliers

teste <- ur.df(capital_trabalho_final, type = "drift", lags = 20, selectlags = "AIC")
summary(teste)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.71942 -0.05589 -0.02009  0.02240  0.84050 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.10407    0.03819   2.725  0.00691 ** 
## z.lag.1     -0.18025    0.06521  -2.764  0.00615 ** 
## z.diff.lag1 -0.51120    0.08142  -6.279 1.61e-09 ***
## z.diff.lag2 -0.33937    0.08193  -4.142 4.78e-05 ***
## z.diff.lag3 -0.31427    0.07891  -3.983 9.06e-05 ***
## z.diff.lag4 -0.32914    0.07457  -4.414 1.54e-05 ***
## z.diff.lag5 -0.16080    0.06372  -2.524  0.01227 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.148 on 238 degrees of freedom
## Multiple R-squared:  0.347,  Adjusted R-squared:  0.3305 
## F-statistic: 21.08 on 6 and 238 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -2.7642 3.8382 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.44 -2.87 -2.57
## phi1  6.47  4.61  3.79

# The value of -2.7642 indicates that we fail to reject the hypothesis of a unit root at 5% sig. level
# Let's take the first difference of this series
capital_trabalho_diff <- diff(capital_trabalho_final)
teste <- ur.df(capital_trabalho_diff, type = "drift", lags = 20, selectlags = "AIC")

# The test statistic of -10,74 makes clear that we reject the null hypothesis of unit root
summary(teste)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73363 -0.05954 -0.01229  0.03755  0.87206 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.001866   0.009584   0.195  0.84579    
## z.lag.1     -3.443603   0.320516 -10.744  < 2e-16 ***
## z.diff.lag1  1.766775   0.284779   6.204 2.43e-09 ***
## z.diff.lag2  1.270467   0.233790   5.434 1.36e-07 ***
## z.diff.lag3  0.822876   0.180244   4.565 8.01e-06 ***
## z.diff.lag4  0.379660   0.124186   3.057  0.00249 ** 
## z.diff.lag5  0.112135   0.064543   1.737  0.08362 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1497 on 237 degrees of freedom
## Multiple R-squared:  0.7665, Adjusted R-squared:  0.7606 
## F-statistic: 129.7 on 6 and 237 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -10.7439 57.7163 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.44 -2.87 -2.57
## phi1  6.47  4.61  3.79

# Compara a série ajustada com a série original e a série diferenciada
df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho_final[1:(length(capital_trabalho_final)-1)]), by = "1 month"), capital_trabalho[1:(length(capital_trabalho)-1)], capital_trabalho_final[1:(length(capital_trabalho_final)-1)], capital_trabalho_diff)
names(df) <- c("Data", "Original", "Série Filtrada", "Série diferenciada")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = capital_trabalho[1:(length(capital_trabalho)-1)], colour = "Original")) +
  geom_line(aes(y = capital_trabalho_final[1:(length(capital_trabalho_final)-1)], colour = "Filtrada")) + 
  geom_line(aes(y = capital_trabalho_diff, colour = "Diferenciada")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
  scale_y_continuous(name = "Razão Capital/Trabalho")+
  labs(color="Variável") +
  theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7))
p1

rm(capital_trabalho_a, capital_trabalho_b, capital_trabalho_2a, capital_trabalho_2b, p1, capital, trabalho, teste)

2. Selic (interest rate)

meses <- ggmonthplot(selic) + 
  theme_bw() +
  scale_y_continuous(name = "Selic (%a.a.)") +
  #labs(title = "Média diária e mensal da taxa Selic", x = "Mês", subtitle = "Nenhuma transformação")+
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))

# Gráfico por ano e mês
season <- ggseasonplot(selic, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Selic (%a.a.)") +
  #labs(title = "Média mensal da taxa Selic", x = "Mês", subtitle = "Nenhuma transformação")+
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Season_SELIC.pdf", width = 6, height = 3)
season

#dev.off()

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Month_SELIC.pdf", width = 6, height = 3)
meses

#dev.off()

bp_ts <- breakpoints(selic ~ 1)
ci_ts <- confint(bp_ts)

Data <- seq(as.Date("1996-01-01"), length = length(selic), by = "1 month")
df1 <- melt(data.frame(Data,selic), id.vars = "Data")
names(df1) <- c("Data", "Variavel", "Valor")

p1 <- ggplot(df1[which(df1$Variavel == "selic"),], aes(Data, Valor, colour = Variavel)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="Selic (%a.a.)") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+ 
        geom_vline(xintercept = Data[bp_ts$breakpoints], linetype="longdash")+
        geom_segment(aes(x = Data[ci_ts$confint[1]], y = min(selic), xend = Data[ci_ts$confint[5]], yend = min(selic)))+
        geom_segment(aes(x = Data[ci_ts$confint[2]], y = min(selic), xend = Data[ci_ts$confint[6]], yend = min(selic)))+
        theme_bw()
p1 <- p1 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), legend.position = "none")

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Quebra_SELIC.pdf", width = 6, height = 3)
p1

#dev.off()

###########
# Dividindo a selic em três pedaços
##########

# Salvando três séries para passar o seasonal
selic_a <- selic[1:(bp_ts$breakpoints[1]-1)]
selic_b <- selic[(bp_ts$breakpoints[1]):(bp_ts$breakpoints[2]-1)]
selic_c <- selic[(bp_ts$breakpoints[2]):(length(selic))]
selic_a <- ts(selic_a,  start = c(1996, 01, 01), frequency = 12)
selic_b <- ts(selic_b,  start = c(1999, 05, 01), frequency = 12)
selic_c <- ts(selic_c,  start = c(2006, 08, 01), frequency = 12)

selic_2a <- selic_a

m <- seas(x = selic_b, transform.function = "none", regression.aictest = NULL)
qs(m)

##                   qs   p-val
## qsori       30.28836 0.00000
## qsorievadj  30.28836 0.00000
## qsrsd        0.00000 1.00000
## qssadj       3.00726 0.22232
## qssadjevadj  3.00726 0.22232
## qsirr        0.00000 1.00000
## qsirrevadj   0.00000 1.00000

summary(m)

## 
## Call:
## seas(x = selic_b, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AR-Nonseasonal-01 -0.19030    0.14589  -1.304    0.192    
## AR-Nonseasonal-02  0.09671    0.09194   1.052    0.293    
## AR-Nonseasonal-03  0.51926    0.08351   6.218 5.03e-10 ***
## MA-Nonseasonal-01  0.12082    0.17475   0.691    0.489    
## MA-Seasonal-12     0.99913    0.12930   7.727 1.10e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (3 1 1)(0 1 1)  Obs.: 87  Transform: none
## AICc: 306.6, BIC: 319.2  QS (no seasonality in final):3.007  
## Box-Ljung (no autocorr.): 44.68 ** Shapiro (normality): 0.9865  
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
##   been found in one or more of the estimated spectra.

selic_2b <- final(m)

m <- seas(x = selic_c, transform.function = "none", regression.aictest = NULL) 
qs(m)

##                    qs   p-val
## qsori        45.47230 0.00000
## qsorievadj   45.47230 0.00000
## qsrsd         0.00000 1.00000
## qssadj        2.90012 0.23456
## qssadjevadj   2.90012 0.23456
## qsirr         0.00000 1.00000
## qsirrevadj    0.00000 1.00000
## qssori       31.25379 0.00000
## qssorievadj  31.25379 0.00000
## qssrsd        0.00000 1.00000
## qsssadj       2.79444 0.24728
## qsssadjevadj  2.79444 0.24728
## qssirr        0.00000 1.00000
## qssirrevadj   0.00000 1.00000

summary(m)

## 
## Call:
## seas(x = selic_c, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AR-Nonseasonal-01 -0.17009    0.10795  -1.576  0.11511    
## AR-Nonseasonal-02  0.21991    0.07042   3.123  0.00179 ** 
## AR-Nonseasonal-03  0.58243    0.06440   9.044  < 2e-16 ***
## MA-Nonseasonal-01  0.29636    0.13073   2.267  0.02339 *  
## MA-Seasonal-12     0.99998    0.09303  10.749  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (3 1 1)(0 1 1)  Obs.: 139  Transform: none
## AICc: 311.3, BIC: 327.6  QS (no seasonality in final):  2.9  
## Box-Ljung (no autocorr.): 61.92 *** Shapiro (normality): 0.9843  
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
##   been found in one or more of the estimated spectra.

selic_2c <- final(m)

###############
# Juntando novamente as três séries
###############

comb <- ts.union(selic_2a, selic_2b, selic_2c)
selic_final <- pmin(comb[,1], comb[,2], comb[,3], na.rm = TRUE)

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(selic_final), by = "1 month"), selic, selic_final)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = selic, colour = "Original"), linetype = "dashed") +
  geom_line(aes(y = selic_final, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Ano", labels = date_format("%Y"))+
  scale_y_continuous(name = "Selic (%a.a.)")+
  labs(color="Série") +
  scale_colour_manual("Variável", breaks = c("Original", "Filtrada"), values = c(cores[1], cores[2]))+
  theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7), legend.position = "top")
p1

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-AntesDepois_SELIC.pdf", width = 6, height = 3)
p1

#dev.off()

teste <- ur.df(selic_final, type = "drift", selectlags = "AIC", lags = 20)
summary(teste)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.4155  -0.8920  -0.1325   0.7859  21.8141 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.55727    0.45884   1.215   0.2258    
## z.lag.1     -0.04258    0.02739  -1.555   0.1213    
## z.diff.lag1 -0.07486    0.06754  -1.108   0.2688    
## z.diff.lag2 -0.12590    0.06436  -1.956   0.0516 .  
## z.diff.lag3 -0.09387    0.06487  -1.447   0.1492    
## z.diff.lag4 -0.09785    0.06484  -1.509   0.1326    
## z.diff.lag5 -0.05325    0.06442  -0.827   0.4092    
## z.diff.lag6  0.05328    0.06387   0.834   0.4051    
## z.diff.lag7 -0.27783    0.06296  -4.413 1.55e-05 ***
## z.diff.lag8 -0.09719    0.06503  -1.495   0.1364    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.594 on 235 degrees of freedom
## Multiple R-squared:  0.1435, Adjusted R-squared:  0.1107 
## F-statistic: 4.376 on 9 and 235 DF,  p-value: 2.692e-05
## 
## 
## Value of test-statistic is: -1.5549 1.4189 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.44 -2.87 -2.57
## phi1  6.47  4.61  3.79

# Eu não acho que faz muito sentido diferenciar a selic pra tirar a raiz unitária... mas é uma possibilidade pra ver se são as variações na selic que impactam nas variações da razão capital trabalho

rm(p1, p2, selic_4390, selic_2a, selic_2b, selic_2c, selic_a, selic_b, selic_c, teste)

2. Swap DI pré 3 meses

meses <- ggmonthplot(swap90) + 
  theme_bw() +
  scale_y_continuous(name = "Swap DI 90d (%a.a)") +
  #labs(title = "Média diária e mensal da taxa Selic", x = "Mês", subtitle = "Nenhuma transformação")+
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))

# Gráfico por ano e mês
season <- ggseasonplot(swap90, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Swap DI 90d (%a.a.)") +
  #labs(title = "Média mensal da taxa swapDI 90d", x = "Mês", subtitle = "Nenhuma transformação")+
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Season_SWAP.pdf", width = 6, height = 3)
season

#dev.off()

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Month_SWAP.pdf", width = 6, height = 3)
meses

#dev.off()

bp_ts <- breakpoints(swap90 ~ 1)
ci_ts <- confint(bp_ts)

#m <- seas(swap90, transform.function = "none", regression.aictest = NULL)
#summary(m) # Series should not be a candidate for seasonal adjustment because the spectrum of the prior adjusted series (Table B1) has no visually significant seasonal peaks.

Data <- seq(as.Date("1999-09-01"), length = length(swap90), by = "1 month")
df1 <- melt(data.frame(Data,swap90), id.vars = "Data")
names(df1) <- c("Data", "Variavel", "Valor")

p1 <- ggplot(df1[which(df1$Variavel == "swap90"),], aes(Data, Valor, colour = Variavel)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="Swap DI 90d (%a.a.)") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+ 
        geom_vline(xintercept = Data[bp_ts$breakpoints], linetype="longdash")+
        geom_segment(aes(x = Data[ci_ts$confint[1]], y = min(swap90), xend = Data[ci_ts$confint[5]], yend = min(swap90)))+
        geom_segment(aes(x = Data[ci_ts$confint[2]], y = min(swap90), xend = Data[ci_ts$confint[6]], yend = min(swap90)))+
        theme_bw()
p1 <- p1 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), legend.position = "none")
p1

## Cortar a selic a partir de 1999-09-01
selic_cortada <- window(selic_final, c(1999, 09))

df8 <- data.frame(Data, swap90, selic_cortada)
names(df8) <- c("Data", "Swap", "Selic")

df9 <- melt(data = df8, id.vars = "Data")

p9 <- ggplot(df9, aes(Data, value, colour = variable))+
  geom_line(alpha = 1, aes(linetype = variable))+
  labs(title="", y = "Taxa (%a.a.)", x = "Ano", color = "Taxa") +
  scale_colour_brewer(palette = "Set1")+
  scale_x_date(date_breaks = "1 year", name = "Ano", labels = date_format("%Y"))+
  theme_bw()

p9 <- p9 + labs(linetype = "Taxa")
p9 <- p9 + labs(colour = "Taxa")
p9 <- p9 + theme(legend.position = "top", legend.key.size = unit(.5, "cm"), axis.text.x = element_text(angle = 25, hjust = 1, size = 7), axis.title.y = element_text(size = 7), axis.title.x = element_text(size = 7))

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Swap_Selic.pdf", width = 6, height = 3)
p9

#dev.off()

3. PIB mensal

#######################
# PIB Mensal          #
#######################
cores <- brewer.pal(6, "Dark2")

meses <- ggmonthplot(PIBmensal) + 
  theme_bw() +
  scale_y_continuous(name = "PIB mensal (R$ - milhões)") +
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))
  #labs(title = "Média diária e mensal da variação do PIB", x = "Mês", subtitle = "Nenhuma transformação")

# Gráfico por ano e mês
season <- ggseasonplot(PIBmensal, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "PIB mensal (R$ - milhões)") +
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))
  #labs(title = "Média mensal da da variação do PIB", x = "Mês", subtitle = "Nenhuma transformação")

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Season_PIB.pdf", width = 6, height = 3)
season

#dev.off()

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Month_PIB.pdf", width = 6, height = 3)
meses

#dev.off()

bp_ts <- breakpoints(PIBmensal ~ 1)
ci_ts <- confint(bp_ts)

Data <- seq(as.Date("1996-01-01"), length = length(PIBmensal), by = "1 month")
df1 <- melt(data.frame(Data,PIBmensal), id.vars = "Data")
names(df1) <- c("Data", "Variavel", "Valor")

p1 <- ggplot(df1[which(df1$Variavel == "PIBmensal"),], aes(Data, Valor, colour = Variavel)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="PIB mensal (R$ milhões)") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+ 
        geom_vline(xintercept = Data[bp_ts$breakpoints], linetype="longdash")+
        geom_segment(aes(x = Data[ci_ts$confint[1]], y = min(PIBmensal), xend = Data[ci_ts$confint[5]], yend = min(PIBmensal)))+
        geom_segment(aes(x = Data[ci_ts$confint[2]], y = min(PIBmensal), xend = Data[ci_ts$confint[6]], yend = min(PIBmensal)))+
        geom_segment(aes(x = Data[ci_ts$confint[3]], y = min(PIBmensal), xend = Data[ci_ts$confint[7]], yend = min(PIBmensal)))+
        geom_segment(aes(x = Data[ci_ts$confint[4]], y = min(PIBmensal), xend = Data[ci_ts$confint[8]], yend = min(PIBmensal)))+
        theme_bw()
p1 <- p1 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), legend.position = "none")

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Quebra_PIB.pdf", width = 6, height = 3)
p1

#dev.off()

m <- seas(x = PIBmensal, transform.function = "none", regression.aictest = NULL)
qs(m)

##                     qs   p-val
## qsori        223.05943 0.00000
## qsorievadj   223.05943 0.00000
## qsrsd          2.04455 0.35978
## qssadj         0.11717 0.94310
## qssadjevadj    0.11717 0.94310
## qsirr          0.00000 1.00000
## qsirrevadj     0.00000 1.00000
## qssori        73.70678 0.00000
## qssorievadj   73.70678 0.00000
## qssrsd         0.65785 0.71970
## qsssadj        0.01014 0.99494
## qsssadjevadj   0.01014 0.99494
## qssirr         0.00000 1.00000
## qssirrevadj    0.00000 1.00000

summary(m)

## 
## Call:
## seas(x = PIBmensal, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AR-Nonseasonal-01 -0.83700    0.13651  -6.131 8.71e-10 ***
## AR-Nonseasonal-02 -0.37993    0.05744  -6.615 3.72e-11 ***
## MA-Nonseasonal-01 -0.57124    0.14263  -4.005 6.20e-05 ***
## MA-Seasonal-12     0.72433    0.04796  15.104  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (2 1 1)(0 1 1)  Obs.: 267  Transform: none
## AICc:  5370, BIC:  5388  QS (no seasonality in final):0.1172  
## Box-Ljung (no autocorr.): 59.13 *** Shapiro (normality): 0.997  
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
##   been found in one or more of the estimated spectra.

# Even though the series presents structural breaks, I'm not dividing to make the adjustment because for the individual pieces the seasonal adjustment is not effective

# Salva a série final em uma nova variável
PIB2 <- final(m)

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(PIBmensal), by = "1 month"), PIBmensal, PIB2)
names(df) <- c("Data", "Original", "Série Filtrada")
p2 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = PIBmensal, colour = "Original"), linetype = "dashed") +
  geom_line(aes(y = PIB2, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Ano", labels = date_format("%Y"))+
  scale_y_continuous(name = "PIB Mensal (R$ milhões)")+
  labs(color="Série") +
  scale_colour_manual("Variável", breaks = c("Original", "Filtrada"), values = c(cores[1], cores[2]))+
  theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7), legend.position = "top")

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-AntesDepois_PIB.pdf", width = 6, height = 3)
p2

#dev.off()

# Now we transform to GDP growth rate
PIBm.df <- ur.df(PIB2, type = "trend", selectlags = "AIC", lags = 20)
summary(PIBm.df)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25069.8  -4087.3     88.1   4655.8  15070.3 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   6.906e+03  3.433e+03   2.011  0.04548 * 
## z.lag.1      -2.806e-02  1.661e-02  -1.690  0.09250 . 
## tt            4.090e+01  2.675e+01   1.529  0.12774   
## z.diff.lag1  -2.028e-01  6.625e-02  -3.061  0.00247 **
## z.diff.lag2  -1.240e-01  6.741e-02  -1.840  0.06713 . 
## z.diff.lag3   4.834e-02  6.779e-02   0.713  0.47652   
## z.diff.lag4  -1.466e-01  6.769e-02  -2.166  0.03140 * 
## z.diff.lag5   1.144e-01  6.827e-02   1.675  0.09523 . 
## z.diff.lag6   1.755e-01  6.793e-02   2.584  0.01040 * 
## z.diff.lag7   3.474e-02  6.808e-02   0.510  0.61033   
## z.diff.lag8   4.757e-02  6.785e-02   0.701  0.48400   
## z.diff.lag9   1.069e-01  6.764e-02   1.580  0.11554   
## z.diff.lag10  8.834e-03  6.804e-02   0.130  0.89681   
## z.diff.lag11  8.502e-02  6.788e-02   1.253  0.21164   
## z.diff.lag12 -2.062e-02  6.812e-02  -0.303  0.76244   
## z.diff.lag13  9.059e-02  6.810e-02   1.330  0.18483   
## z.diff.lag14  2.891e-02  6.799e-02   0.425  0.67113   
## z.diff.lag15  1.005e-01  6.871e-02   1.463  0.14482   
## z.diff.lag16 -5.558e-02  6.828e-02  -0.814  0.41655   
## z.diff.lag17  1.853e-01  6.837e-02   2.710  0.00724 **
## z.diff.lag18  1.134e-01  6.839e-02   1.659  0.09855 . 
## z.diff.lag19 -1.103e-01  6.716e-02  -1.642  0.10192   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6869 on 224 degrees of freedom
## Multiple R-squared:  0.299,  Adjusted R-squared:  0.2333 
## F-statistic:  4.55 on 21 and 224 DF,  p-value: 2.706e-09
## 
## 
## Value of test-statistic is: -1.6896 1.9696 1.5483 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

# Since it is non-stationary, we take log difference
PIB <- diff(log(PIB2), 1)
autoplot(PIB)

PIB.df <- ur.df(PIB, type = "none", selectlags = "AIC", lags = 20)
summary(PIB.df)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.057193 -0.007917  0.002290  0.010353  0.046704 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## z.lag.1      -0.42364    0.24624  -1.720 0.086730 .  
## z.diff.lag1  -0.75572    0.24622  -3.069 0.002410 ** 
## z.diff.lag2  -0.87318    0.24710  -3.534 0.000498 ***
## z.diff.lag3  -0.92117    0.24993  -3.686 0.000286 ***
## z.diff.lag4  -1.13965    0.25271  -4.510 1.05e-05 ***
## z.diff.lag5  -1.02745    0.25598  -4.014 8.15e-05 ***
## z.diff.lag6  -0.90695    0.25810  -3.514 0.000534 ***
## z.diff.lag7  -0.90280    0.25660  -3.518 0.000526 ***
## z.diff.lag8  -0.83808    0.25418  -3.297 0.001136 ** 
## z.diff.lag9  -0.70885    0.24920  -2.845 0.004859 ** 
## z.diff.lag10 -0.72064    0.24341  -2.961 0.003400 ** 
## z.diff.lag11 -0.65663    0.23631  -2.779 0.005922 ** 
## z.diff.lag12 -0.69323    0.22812  -3.039 0.002657 ** 
## z.diff.lag13 -0.64869    0.21712  -2.988 0.003124 ** 
## z.diff.lag14 -0.63229    0.20346  -3.108 0.002130 ** 
## z.diff.lag15 -0.54885    0.18955  -2.896 0.004160 ** 
## z.diff.lag16 -0.63232    0.17031  -3.713 0.000259 ***
## z.diff.lag17 -0.41204    0.14433  -2.855 0.004711 ** 
## z.diff.lag18 -0.31236    0.11932  -2.618 0.009452 ** 
## z.diff.lag19 -0.34610    0.09137  -3.788 0.000195 ***
## z.diff.lag20 -0.23364    0.06153  -3.797 0.000188 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01582 on 224 degrees of freedom
## Multiple R-squared:  0.7146, Adjusted R-squared:  0.6878 
## F-statistic:  26.7 on 21 and 224 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -1.7205 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

# Let's try to plot everything together
df <- data.frame(seq(as.Date("1996-01-01"), length = length(PIBmensal[-length(PIBmensal)]), by = "1 month"), PIBmensal[-length(PIBmensal)]/10^(6), PIB2[-length(PIB2)]/10^(6), PIB)
names(df) <- c("Data", "Original", "Série Filtrada", "Série diferenciada")
p3 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = PIBmensal[-length(PIBmensal)]/10^(6), colour = "Original")) +
  geom_line(aes(y = PIB2[-length(PIB2)]/10^(6), colour = "Filtrada")) + 
  geom_line(aes(y = PIB, colour = "Diferenciada")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
  scale_y_continuous(name = "PIB Mensal (base = 03/18)")+
  labs(color="Variável") +
  theme_bw()
p3 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 8))

4. Exchange rate

#######################
# Câmbio              #
#######################
cores <- brewer.pal(6, "Dark2")

meses <- ggmonthplot(cambio_raw) + 
  theme_bw() +
  scale_y_continuous(name = "Taxa de Câmbio Efetiva Real (índice)") +
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))
  #labs(title = "Média diária e mensal da taxa de câmbio", x = "Month", subtitle = "No data transformation")

# Gráfico por ano e mês
season <- ggseasonplot(cambio_raw, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Taxa de Câmbio Efetiva Real (índice)") +
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))
  #labs(title = "Monthly average of exchange rate", x = "Month", subtitle = "No data transformation")

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Season_CAMBIO.pdf", width = 6, height = 3)
season

#dev.off()

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Month_CAMBIO.pdf", width = 6, height = 3)
meses

#dev.off() 

bp_ts <- breakpoints(cambio_raw ~ 1)
ci_ts <- confint(bp_ts)

Data <- seq(as.Date("1996-01-01"), length = length(cambio_raw), by = "1 month")
df1 <- melt(data.frame(Data,cambio_raw), id.vars = "Data")
names(df1) <- c("Data", "Variavel", "Valor")

p1 <- ggplot(df1[which(df1$Variavel == "cambio_raw"),], aes(Data, Valor, colour = Variavel)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="Taxa de Câmbio Efetiva Real (índice)") +
        scale_x_date(date_breaks = "12 months", name = "Date", labels = date_format("%m-%Y"))+ 
        geom_vline(xintercept = Data[bp_ts$breakpoints], linetype="longdash")+
        geom_segment(aes(x = Data[ci_ts$confint[1]], y = min(cambio_raw), xend = Data[ci_ts$confint[5]], yend = min(cambio_raw)))+
        geom_segment(aes(x = Data[ci_ts$confint[2]], y = min(cambio_raw), xend = Data[ci_ts$confint[7]], yend = min(cambio_raw)))+
        geom_segment(aes(x = Data[ci_ts$confint[3]], y = min(cambio_raw), xend = Data[ci_ts$confint[8]], yend = min(cambio_raw)))+
        geom_segment(aes(x = Data[ci_ts$confint[4]], y = min(cambio_raw), xend = Data[ci_ts$confint[9]], yend = min(cambio_raw)))+
        geom_segment(aes(x = Data[ci_ts$confint[5]], y = min(cambio_raw), xend = Data[ci_ts$confint[10]], yend = min(cambio_raw)))+
        theme_bw()
p1 <- p1 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), legend.position = "none")
p1

m <- seas(x = cambio_raw, transform.function = "none", regression.aictest = NULL)
qs(m)

##                   qs   p-val
## qsori        0.00000 1.00000
## qsorievadj   0.00000 1.00000
## qsrsd        0.00000 1.00000
## qssadj       0.00000 1.00000
## qssadjevadj  0.00000 1.00000
## qsirr        0.57866 0.74876
## qsirrevadj   0.00000 1.00000
## qssori       0.00000 1.00000
## qssorievadj  0.00000 1.00000
## qssrsd       0.00000 1.00000
## qsssadj      0.00000 1.00000
## qsssadjevadj 0.00000 1.00000
## qssirr       0.00091 0.99954
## qssirrevadj  0.00091 0.99954

summary(m)

## 
## Call:
## seas(x = cambio_raw, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## AO1999.Jan         17.02417    3.10157   5.489 4.04e-08 ***
## LS1999.Feb         38.41276    5.45786   7.038 1.95e-12 ***
## AO2002.Oct         13.48199    1.51561   8.895  < 2e-16 ***
## AO2003.Jan        -10.34927    1.51561  -6.828 8.58e-12 ***
## AO2008.Dec          7.21518    1.47400   4.895 9.83e-07 ***
## MA-Nonseasonal-01  -0.54829    0.05163 -10.619  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (0 1 1)  Obs.: 266  Transform: none
## AICc:  1367, BIC:  1391  QS (no seasonality in final):    0  
## Box-Ljung (no autocorr.): 30.11   Shapiro (normality): 0.9482 ***
## Messages generated by X-13:
## Warnings:
## - At least one visually significant seasonal peak has been
##   found in one or more of the estimated spectra.
## 
## Notes:
## - Unable to test LS1999.Jan due to regression matrix
##   singularity.
## - Unable to test LS1999.Jan due to regression matrix
##   singularity.
## - Unable to test LS2015.Aug due to regression matrix
##   singularity.
## - Unable to test LS2008.Sep due to regression matrix
##   singularity.

cambio <- cambio_raw
# Não tem evidência que precisa ajuste sazonal na série toda
# Como não fiz nenhuma transformação não tem muito o que plotar de gráfico aqui

5. IPCA

#######################
# IPCA                #
#######################
cores <- brewer.pal(6, "Dark2")

meses <- ggmonthplot(ipca) + 
  theme_bw() +
  scale_y_continuous(name = "Taxa de Inflação (IPCA - % a.a.)") +
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))
  #labs(title = "Monthly average of the inflation index", x = "Month", subtitle = "No data transformation")

# Gráfico por ano e mês
season <- ggseasonplot(ipca, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Taxa de Inflação (IPCA - % a.a.)") +
  labs(title = "", x = "Mês", subtitle = "")+
  theme(axis.title = element_text(size=7), axis.text.x = element_text(angle=25, hjust = 1, size = 7))
  #labs(title = "Monthly average of the inflation index", x = "Month", subtitle = "No data transformation")

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Season_IPCA.pdf", width = 6, height = 3)
season

#dev.off()

#pdf(file = "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Imagens\\Fig-Month_IPCA.pdf", width = 6, height = 3)
meses

#dev.off()


bp_ts <- breakpoints(ipca ~ 1)
ci_ts <- confint(bp_ts)

Data <- seq(as.Date("1996-01-01"), length = length(ipca), by = "1 month")
df1 <- melt(data.frame(Data,ipca), id.vars = "Data")
names(df1) <- c("Data", "Variavel", "Valor")

p1 <- ggplot(df1[which(df1$Variavel == "ipca"),], aes(Data, Valor, colour = Variavel)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="Inflation (% monthly)") +
        scale_x_date(date_breaks = "12 months", name = "Date", labels = date_format("%m-%Y"))+ 
        geom_vline(xintercept = Data[bp_ts$breakpoints], linetype="longdash")+
        geom_segment(aes(x = Data[ci_ts$confint[1]], y = min(ipca), xend = Data[ci_ts$confint[5]], yend = min(ipca)))+
        geom_segment(aes(x = Data[ci_ts$confint[2]], y = min(ipca), xend = Data[ci_ts$confint[6]], yend = min(ipca)))+
        theme_bw()
p1 <- p1 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), legend.position = "none")
p1

# Now I divide the series in 3 different pieces to apply the seasonal filter
# Salvando três séries para passar o seasonal
ipca_a <- ipca[1:(bp_ts$breakpoints[1]-1)]
ipca_b <- ipca[(bp_ts$breakpoints[1]):(bp_ts$breakpoints[2]-1)]
ipca_c <- ipca[(bp_ts$breakpoints[2]):(length(ipca))]
ipca_a <- ts(ipca_a,  start = c(1996, 01, 01), frequency = 12)
ipca_b <- ts(ipca_b,  start = c(2002, 03, 01), frequency = 12)
ipca_c <- ts(ipca_c,  start = c(2005, 06, 01), frequency = 12)

# There is no evidence of seasonal component in either the complete series nor the pieces
m <- seas(x = ipca_a, transform.function = "none", regression.aictest = NULL)
qs(m)

##                  qs   p-val
## qsori       0.00000 1.00000
## qsorievadj  0.00000 1.00000
## qsrsd       0.52144 0.77050
## qssadj      0.00000 1.00000
## qssadjevadj 0.00000 1.00000
## qsirr       0.00000 1.00000
## qsirrevadj  0.00000 1.00000
## qssrsd      0.35015 0.83939

summary(m)

## 
## Call:
## seas(x = ipca_a, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AR-Nonseasonal-01  0.79459    0.06821  11.649   <2e-16 ***
## MA-Seasonal-12     0.50638    0.09393   5.391    7e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (1 1 0)(0 0 1)  Obs.: 74  Transform: none
## AICc: 92.27, BIC: 98.79  QS (no seasonality in final):    0  
## Box-Ljung (no autocorr.): 22.07   Shapiro (normality): 0.984  
## Messages generated by X-13:
## Warnings:
## - Series should not be a candidate for seasonal adjustment
##   because the spectrum of the original series (Table A1 or
##   B1) has no visually significant seasonal peaks.
## 
## Notes:
## - Model used for SEATS decomposition is different from the
##   model estimated in the regARIMA modeling module of
##   X-13ARIMA-SEATS.

# Não tem evidência que precisa ajuste sazonal na primeira parte
m <- seas(x = ipca_b, transform.function = "none", regression.aictest = NULL)
qs(m)

##             qs p-val
## qsori        0     1
## qsorievadj   0     1
## qsrsd        0     1
## qssadj       0     1
## qssadjevadj  0     1
## qsirr        0     1
## qsirrevadj   0     1

summary(m)

## 
## Call:
## seas(x = ipca_b, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## LS2002.Nov          1.4617     0.3866   3.781 0.000156 ***
## LS2003.Nov         -1.5678     0.3866  -4.056    5e-05 ***
## AR-Nonseasonal-01   0.7847     0.0954   8.226  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (1 1 0)  Obs.: 39  Transform: none
## AICc: 64.01, BIC: 69.35  QS (no seasonality in final):    0  
## Box-Ljung (no autocorr.):  33.2 . Shapiro (normality): 0.9852  
## Messages generated by X-13:
## Warnings:
## - Series should not be a candidate for seasonal adjustment
##   because the spectrum of the prior adjusted series (Table
##   B1) has no visually significant seasonal peaks.
## 
## Notes:
## - The program cannot perform hypothesis tests for kurtosis
##   on less than 50 observations.

# Não tem evidência que precisa ajuste sazonal na segunda parte
m <- seas(x = ipca_c, transform.function = "none", regression.aictest = NULL)
qs(m)

##                   qs   p-val
## qsori        0.00000 1.00000
## qsorievadj   0.00000 1.00000
## qsrsd        2.06009 0.35699
## qssadj       0.00000 1.00000
## qssadjevadj  0.00000 1.00000
## qsirr        0.00000 1.00000
## qsirrevadj   0.00000 1.00000
## qssori       0.00000 1.00000
## qssorievadj  0.00000 1.00000
## qssrsd       1.09881 0.57729
## qsssadj      0.00000 1.00000
## qsssadjevadj 0.00000 1.00000
## qssirr       0.00000 1.00000
## qssirrevadj  0.00000 1.00000

summary(m)

## 
## Call:
## seas(x = ipca_c, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AR-Nonseasonal-01  0.63263    0.05707   11.09   <2e-16 ***
## MA-Seasonal-12     0.99989    0.06319   15.82   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (1 1 0)(0 0 1)  Obs.: 153  Transform: none
## AICc: -19.97, BIC: -11.06  QS (no seasonality in final):    0  
## Box-Ljung (no autocorr.): 34.49 . Shapiro (normality): 0.993  
## Messages generated by X-13:
## Warnings:
## - Series should not be a candidate for seasonal adjustment
##   because the spectrum of the original series (Table A1 or
##   B1) has no visually significant seasonal peaks.
## 
## Notes:
## - Model used for SEATS decomposition is different from the
##   model estimated in the regARIMA modeling module of
##   X-13ARIMA-SEATS.

# Não tem evidência que precisa ajuste sazonal na terceira parte

m <- seas(x = ipca, transform.function = "none", regression.aictest = NULL)
qs(m)

##                   qs   p-val
## qsori        0.00000 1.00000
## qsorievadj   0.00000 1.00000
## qsrsd        1.84107 0.39830
## qssadj       0.00000 1.00000
## qssadjevadj  0.00000 1.00000
## qsirr        0.00000 1.00000
## qsirrevadj   0.00000 1.00000
## qssori       0.00000 1.00000
## qssorievadj  0.00000 1.00000
## qssrsd       2.06590 0.35595
## qsssadj      0.00000 1.00000
## qsssadjevadj 0.00000 1.00000
## qssirr       0.00000 1.00000
## qssirrevadj  0.00000 1.00000

summary(m)

## 
## Call:
## seas(x = ipca, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## LS2003.Jul        -0.55724    0.12823  -4.346 1.39e-05 ***
## AR-Nonseasonal-01  0.76739    0.03792  20.239  < 2e-16 ***
## MA-Seasonal-12     0.89032    0.02928  30.404  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (1 1 0)(0 0 1)  Obs.: 266  Transform: none
## AICc:   179, BIC: 193.2  QS (no seasonality in final):    0  
## Box-Ljung (no autocorr.):    14   Shapiro (normality): 0.9356 ***
## Messages generated by X-13:
## Warnings:
## - Series should not be a candidate for seasonal adjustment
##   because the spectrum of the prior adjusted series (Table
##   B1) has no visually significant seasonal peaks.
## 
## Notes:
## - Model used for SEATS decomposition is different from the
##   model estimated in the regARIMA modeling module of
##   X-13ARIMA-SEATS.

# Não tem evidência que precisa ajuste sazonal na série toda

# Como não fiz nenhuma transformação não tem muito o que plotar de gráfico aqui

Variávels finais após ajuste sazonal

Descritivas

	Capital_trabalho	Swap	PIB	Cambio	IPCA
Observações	266.0000	266.0000	266.0000	266.0000	266.0000
Mínimo	-8.1395	6.5400	-6.5307	66.5600	1.6456
1o quartil	44.8218	10.7125	-0.9087	81.0700	4.8092
Média	55.7700	13.9905	0.2945	97.4647	6.8280
Mediana	51.2403	13.1600	0.3725	95.6100	6.2707
3o quartil	63.6885	17.3900	1.6253	108.0600	7.6704
Máximo	151.0014	27.4900	5.3142	171.0300	21.9876
Desv. Pad.	19.0722	4.7089	1.8553	21.6786	3.4314

Estimating the BVAR

# Junta o swap com a selic
selic_trim <- ts(selic_final[1:44], start = c(1996, 01, 01), frequency = 12)

# Binding the series
comb <- ts.union(selic_trim, swap90)
swap <- pmin(comb[,1], comb[,2], na.rm = TRUE)

#Inflação - Produto - Juros - Câmbio - Desigualdade.
var1 <- cbind(capital_trabalho_final*100, ipca, PIB*100, swap, cambio) # The bvar function does not allows data.frames
#var2 <- cbind(capital_trabalho_final, swap, PIB2, cambio, ipca) # The bvar function does not allows data.frames
names(var1) <- c("capital_trabalho", "ipca", "pib", "swap", "cambio")
#names(var1) <- c("ipca", "pib", "swap", "cambio","capital_trabalho") # Ficou uma merda

#var1 <- c(1,2,3,4,5)

pos_ct     <- grep("capital_trabalho", names(var1))
pos_swap   <- grep("swap", names(var1))
pos_pib    <- grep("pib", names(var1))
pos_cambio <- grep("cambio", names(var1))
pos_ipca   <- grep("ipca", names(var1))


set.seed(6969)

nburn.    <- 1000000 # Burn the first 100kk replications
nrep.     <- 100000 # Estimate another 100k times after burn in
thinfac.  <- 4 # But keep only 1 in 4 estimates
p.        <- 1 # Number of lags
tau.      <- 48 - p. # Number of observations to use as prior
nf.       <- 1 # Number of forecast periods
k_B.      <- 6 # Fator de inflação na priori
k_A.      <- 4 # Fator de inflação na priori

# No idea what's this for
tamanho <- floor(1/thinfac.*nrep.)
# número de tempos
tempo <- nrow(var1)

pm <- proc.time()
#fit1 <- bvar.sv.tvp(var1, p=p., tau = tau., nburn = nburn., nrep = nrep., thinfac = thinfac., nf = nf., k_B = k_B., k_A = k_A.)
#saveRDS(fit1, "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Testes\\fit.rds") 
fit1 <- readRDS("D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Testes\\Usados no artigo\\fit.rds") 
#rm(fit1, ira)
proc.time() - pm

##    user  system elapsed 
##   17.81    1.23   19.08

Caso queira ver qual a priori, descomente esses blocos de código (que estão ocultos).

Funções impulso resposta

Peguei a explicação daqui daqui: A FIR estima o impacto de um choque unitário em algum elemento de $\varepsilon_t$. Os elementos de $c_t$, $B_{j,t}$ e $A_t$ em (1) são mantidos fixos em seus valores em $t$, mas os elementos de $\Sigma_t$ são substituídos pela média de todos os períodos. Note que a ordem das variáveis importa e deve ser justificada por argumentos econômicos. As FIR contém os percentis 5, 25, 50, 75 e 95.

Impactos dos choques monetários

Em todas variáveis do modelo, no último período

Na razão capital trabalho

Na razão capital trabalho porém em períodos de mudança no Banco Central.

Os presidentes do BCB no período analisado foram:

Arminio Fraga Neto - 4/3/1999 a 1/1/2003
Henrique de Campos Meirelles - 1/1/2003 a 1/1/2011 $t=37$
Alexandre Antonio Tombini - 1/1/2011 a 8/6/2016 $t=133$
Ilan Goldfajn - 9/6/2016 - atual $t=198$

Supondo que existe algum tipo de coordenação entre política fiscal e monetária, peguei os Ministros da Fazenda no Período também:

Pedro Malan - 01/01/1995 a 01/01/2003
Antonio Palocci Filho - 01/01/2003 a 27/03/2006 $t=37$
Guido Mantega - 27/03/2006 a 01/01/2015 $t=75$
Joaquim Levy - 01/01/2015 a 18/12/2015 $t=183$
Nelson Barbosa - 18/12/2015 a 12/05/2016 $t=192$
Henrique Meirelles - 12/05/2016 a 10/04/2018 $t=197$
Eduardo Guardia - 10/04/2018 - atual

No PIB

No PIB, de acordo com mudanças no Banco Central ou ministério da Fazenda.

No Câmbio

Impacto dos choques monetários no Câmbio, de acordo com as mudanças no Banco Central ou Ministério da Fazenda.

No IPCA

Impacto dos choques monetários no IPCA, de acordo com as mudanças no Banco Central ou Ministério da Fazenda.

Choques do PIB

Na razão capital trabalho

No Câmbio

No swap

No IPCA

Choques na taxa de câmbio

Efeito do câmbio na razão capital-trabalho

Efeito do câmbio No Swap

Efeito do câmbio No IPCA

Efeito do câmbio No PIB

Choques no IPCA

Efeito do Ipca na razão capital-trabalho

Efeito do Ipca no swap

Efeito do IPCA na taxa de câmbio

Efeito do IPCA no PIB

Avaliando os parâmetros

Os gráficos abaixo contém os betas (esq) e volatilidade (dir) ao longo do tempo, com o 16o e o 68o percentis (como foi feito no artigo do Primiceri).

Razão Capital-Trabalho

Coeficiente (da própria variável defasada) e desvio padrão do resíduo, usando a média e os quantis 16 e 84 (os mesmos do primiceri) - esquerda - e mesmo gráfico, porém com o intervalo entre o 5 e 95 (área azul) quantis e 25 e 75 quantis (área cinza):

Swap

Coeficiente (da própria variável defasada) e desvio padrão do resíduo, usando a média e os quantis 16 e 84 (os mesmos do primiceri) - primeira coluna - e mesmo gráfico, porém com o intervalo entre o 5 e 95 (área azul) quantis e 25 e 75 quantis (área cinza) - segunda coluna. Também coloquei um gráfico do coeficiente com a escala vertical incluindo o zero - terceira coluna - para dar comparabilidade:

PIB

Taxa de Câmbio

IPCA

Plotando as variáveis todas juntas

Primeiro, plotando os gráficos com o intervalo do 16o ao 84o quantil (como feito no Primiceri e como eu fiz no trabalho de Macro2:

Mas o que acontece se incluímos o zero no eixo vertical dos coeficientes?

Agora vamos repetir gráfico mas fazendo as bandas do quantil 5 a 95 (em azul) e 25 a 75 (em cinza) - sem forçar o zero no eixo dos coeficientes:

Gráfico dos coeficientes das outras variáveis defasadas na equação da razão capital-trabalho

Primeiro os gráficos com os intervalos de 16 a 84% sem forçar o zero no eixo y:

Agora o mesmo gráfico mas forçando o zero a aparecer:

E o mesmo gráfico mas com os quantis de 5 a 95 e 25 a 75:

Repete o mesmo de antes, mas incluindo zero no eixo y:

Diagnóstico

Estou me baseando nos seguintes linkis: * www.patricklam.org/teaching/convergence_print.pdf para traceplot, density plot, running mean plot (é o plot da iteração contra a média das observações até a iteração (ele está indo de 10 em 10 para frente na média)) * https://theoreticalecology.wordpress.com/2011/12/09/mcmc-chain-analysis-and-convergence-diagnostics-with-coda-in-r/

Diagnósticos para as médias

Um dos diagnósticos é plotar a cadeia (usando plot(mcmc_object)). Ele retorna o gráfico dos valores (um gráfico de linha - quanto mais ruidoso melhor) e um traceplot que é o gráfico da distribuição. Se eu não estou enganada, os meus Betas tem distribuição à posteriori normal e falta ver a distribuição dos sigmas. Se a densidade das matrizes é inversa wishart eles tem distribuição gamma, mas acho que não tem fórmula fechada para $\Sigma$ - preciso ver no artigo do Primeri depois.

Outro diagnóstico das médias envolve fazer o gráfico da média usando uma janela móvel para ver se a posteriori cai em alguma região e fica empacada.

Como eu tenho muitos instantes de tempo (e portanto uma quantidade enorme de betas e de sigmas), eu selecionei arbitrariamente 15 betas (que é a divisão inteira de t por 10 - então estou plotando 10% dos parâmetros para cada $y$). Esses gráficos eu resolvi deixar de fora do artigo.

## [1] 190

## [1] 190

## [1] 190

## [1] 190

## [1] 190

## [1] 190

## [1] 190

## [1] 190

## [1] 190

## [1] 190

Agora tem o cumuplot() que faz o gráfico dos quantis da cadeia (de forma acumulada). O default é 0.025, 0.5, 0.975. Se eu tiver tempo depois, vou mexer no layout desses gráficos pois são feiosos demais. Também deixei de fora.

Por fim, tem o diagnóstico do Geweke, que calcula uma estatística com distribuição normal e que compara o início e o final da distribuição. De novo vou plotar só uma parte para não ficar muito grande. E deixei de fora do artigo

Diagnóstico da autocorrelação

O bloco abaixo eu montei para entender melhor o que a autocorr() do pacote coda faz. No fim, vou ficar com a função de autocorrelação (não a parcial) pois ela avalia o efeito até o instante j. (o bloco agora está oculto)

Agora estou calculando a função te autocorrelação para depois fazer uns fancy graphs. O bloco do código está oculto.

##    user  system elapsed 
##  109.08    0.01  109.22

Agora vamos para os gráficos…. a ideia é fazer algo como o Primiceri tem no anexo B (figura 9 - topo), porém dividindo por variável do VAR.

Diagnóstico multivariado

A estatística de Gelman-Rubin mede se há diferença significativa entre a variância dentro e entre diversas cadeias utilizando um fator de redução de escala. Aqui eu optei por simular duas cadeias e fazer a comparação entre elas. Em função da memória do computador eu utilizei um número menor de iterações salvas: ao invés de salvar todos os valores, armazenei apenas os 10.000 últimos valores (20%).

##    user  system elapsed 
##   17.44    1.06   20.62

## [1] 1.495861

##   parametro variable    value
## 1         1   gelman 1.001433
## 2         2   gelman 1.001451
## 3         3   gelman 1.001453
## 4         4   gelman 1.001473
## 5         5   gelman 1.001477
## 6         6   gelman 1.001587

Previsão

Fiz previsão $nf = 1, 6, 12, 24$ passos à frente.

O funcionamento do algoritmo the previsão é o seguinte:

Para cada uma das iterações MCMC (no meu caso, $\ell = 1, \ldots, 100.000$);
1. Simula a trajetória dos parâmetros para $nf$ períodos à frente, baseado na iteração atual $\ell$;
2. Simula os vetores de erros $\varepsilon_{T+j}$, $j=1, \ldots, nf$ a partir de uma normal;
3. Usa os parâmetros e erros para simular observações previstas usando o processo gerador dos dados que está no início do documento;
4. Para cada período $j = 1, \ldots, nf$ calcula o vetor de médias e a matriz de variâncias e covariâncias condicional aos parâmetros simulados até $T+j$.
“Encolhe as observações”: fica com uma em cada thinfac observações (aqui usei uma em 4);
Entrega os resultados:
- fc.ydraws é um objeto de dimensão [5 \times nf \times 25.000] (5 variáveis no var, nf é o número de períodos previsto e 25.000 são o número de observações que eu guardo das 100.000 simulações). Por exemplo, fc.ydraws[1,3,] contém todas as repetições da primeira variável do VAR, 3 períodos à frente;
- fc.mdraws contém as médias (que que é isso?) condicionais aos parâmetrs. Seu tamanho é o mesmo de fc.ydraws;
- fc.vdraws contém as variâncias previstas, $V_{T+j}^{(\ell)}$ e tem tamanho $[0.5\cdot 5 \cdot (5+1) \times nf \times 25000]$. Por exemplo, fc.vdraws[,1,1] contém os elementos empilhados de $V_{T+1}^{(1)}$ e é a primeira amostra das matriz de variância e covariância para todas variáveis, um período à frente.

Medidas de erro de previsão usadas

A primeira coisa é implementar uma função para calcular o MASE (mean absolute scaled error). Defina o termo de erro de um passo à frente como $e_t = y_t - \hat{y}_t$. Então, podemos definir um erro padronizado como

\[q_t = \frac{e_t}{\frac{1}{n-1}\sum\limits_{i=2}^n |y_i - y_{i-1}|},\] que é independente da escala dos dados das séries. Então, o erro absoluto padronizado médio (algumas pessoas traduzem como escalado ou escalar, eu achei feio demais) é definido por

\[\text{MASE} = \text{média }(|q_t|)\]

Além desse, vou calcular também outras 2 medidas comuns, MAPE (erro médio percentual absoluto) e o MSE (erro quadrático médio). Esses estão implementados no R dentro do pacote DescTools.

Então vou fazer umas funções auxiliares para calcular o MASE e também para juntar os resultados em uma tabelinha.

# O denominador do MASE é fixo, então vou calcular ele separadamente
# Como eu junto todas as variáveis do VAR no objeto var1, vou aplicar a função sobre as colunas do var1

den_MASE <- function(serie){
  n <- length(serie)
  soma <- vector()
  for (i in 2:n){
    soma[i-1] <- abs(serie[i] - serie[i-1])
  }
  x <- sum(soma)/(n-1)
  x
}

1 passo à frente

~Posso usar o que eu já tinha estimado.~ Não posso usar o que eu tinha estimado porque eu não separei o pedaço da amostra que vai ficar de comparação.

rm(fit1)
# Previsão
#Inflação - Produto - Juros - Câmbio - Desigualdade.
var1 <- cbind(capital_trabalho_final*100, ipca, PIB*100, swap, cambio) # The bvar function does not allows data.frames
#var2 <- cbind(capital_trabalho_final, swap, PIB2, cambio, ipca) # The bvar function does not allows data.frames
names(var1) <- c("capital_trabalho", "ipca", "pib", "swap", "cambio")
#names(var1) <- c("ipca", "pib", "swap", "cambio","capital_trabalho") # Ficou uma merda

#var1 <- c(1,2,3,4,5)

pos_ct     <- grep("capital_trabalho", names(var1))
pos_swap   <- grep("swap", names(var1))
pos_pib    <- grep("pib", names(var1))
pos_cambio <- grep("cambio", names(var1))
pos_ipca   <- grep("ipca", names(var1))

nburn.    <- 1000000 # Burn the first 100kk replications
nrep.     <- 100000 # Estimate another 100k times after burn in
thinfac.  <- 4 # But keep only 1 in 4 estimates
p.        <- 1 # Number of lags
tau.      <- 48 - p. # Number of observations to use as prior
nf.       <- 1 # Number of forecast periods
k_B.      <- 6 # Fator de inflação na priori
k_A.      <- 4 # Fator de inflação na priori

# No idea what's this for
tamanho <- floor(1/thinfac.*nrep.)
# número de tempos
tempo <- nrow(var1)

# Prepara um vetor com as últimas nf. observações para deixar de fora
amostra_1 <- var1[(nrow(var1)-nf.+1):(nrow(var1)),] 
var1      <- var1[1:(nrow(var1)-nf.),]

set.seed(2747)
fit1 <- readRDS("D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Testes\\Usados no artigo\\fit1for.rds")
#fit1 <- bvar.sv.tvp(var1, p=p., tau = tau., nburn = nburn., nrep = nrep., nf = nf., thinfac = thinfac.)
#saveRDS(fit1, "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Testes\\fit1for.rds")
#rm(fit1)

# Primeiro eu monto um objeto só com a previsão
# Eu sei que essa coisa tem três dimensões e quero calcular a média da terceira dimensão e ficar com uma matriz
var1_prev <- apply(fit1$fc.ydraws, c(1,2), mean)

# Agora preciso calcular os erros de previsão em cada período
# nf é o número de períodos de previsão
# amostra1 é onde ficam guardados os valores que estou prevendo
# var1_prev é onde ficam guardados os valores que eu previ

amostra1 <- matrix(amostra_1, ncol=nf.)
erro_prev <- amostra1 - var1_prev

q_t <- erro_prev/apply(var1[((tau.+p.+1):nrow(var1)),], 2, FUN = den_MASE)

MASE_1 <- abs(q_t)

MAPE_1 <- vector()
MSE_1  <- vector()

for (i in 1:ncol(var1)){
   MAPE_1[i] <- MAPE(x = var1_prev[i,1], ref = amostra1[i,1])
   MSE_1[i] <- MSE(x = var1_prev[i,1], ref = amostra1[i,1])
}

Agora estimo um modelo VAR(1) simples para comparar

fit_simples <- VAR(var1, p=p.)
var_simples_prev <- predict(fit_simples, n.ahead = nf.)
vars_prev <- vector()

for (i in 1:ncol(var1)){
  vars_prev[i] <-  (var_simples_prev[[1]][[i]][1])
}

erro_simples <- amostra_1 - vars_prev

q_t <- erro_simples/apply(var1[((tau.+p.+1):nrow(var1)),], 2, FUN = den_MASE)

MASE_1_simples <- abs(q_t)

MAPE_1_simples <- vector()
MSE_1_simples  <- vector()
vars_prev <- matrix(vars_prev, ncol = 1)
amostra_1 <- matrix(amostra_1, ncol = 1)

for (i in 1:ncol(var1)){
   MAPE_1_simples[i] <- MAPE(x = vars_prev[i,1], ref = amostra_1[i,1])
   MSE_1_simples[i] <- MSE(x = vars_prev[i,1], ref = amostra_1[i,1])
}

6 passos à frente (1 semestre)

rm(fit1)
# Previsão
#Inflação - Produto - Juros - Câmbio - Desigualdade.
var1 <- cbind(capital_trabalho_final*100, ipca, PIB*100, swap, cambio) # The bvar function does not allows data.frames
#var2 <- cbind(capital_trabalho_final, swap, PIB2, cambio, ipca) # The bvar function does not allows data.frames
names(var1) <- c("capital_trabalho", "ipca", "pib", "swap", "cambio")
#names(var1) <- c("ipca", "pib", "swap", "cambio","capital_trabalho") # Ficou uma merda

#var1 <- c(1,2,3,4,5)

pos_ct     <- grep("capital_trabalho", names(var1))
pos_swap   <- grep("swap", names(var1))
pos_pib    <- grep("pib", names(var1))
pos_cambio <- grep("cambio", names(var1))
pos_ipca   <- grep("ipca", names(var1))

nburn.    <- 1000000 # Burn the first 100kk replications
nrep.     <- 100000 # Estimate another 100k times after burn in
thinfac.  <- 4 # But keep only 1 in 4 estimates
p.        <- 1 # Number of lags
tau.      <- 48 - p. # Number of observations to use as prior
nf.       <- 6 # Number of forecast periods
k_B.      <- 6 # Fator de inflação na priori
k_A.      <- 4 # Fator de inflação na priori

# No idea what's this for
tamanho <- floor(1/thinfac.*nrep.)
# número de tempos
tempo <- nrow(var1)

# Prepara um vetor com as últimas nf. observações para deixar de fora
amostra_1 <- var1[(nrow(var1)-nf.+1):(nrow(var1)),] 
var1      <- var1[1:(nrow(var1)-nf.),]

set.seed(2747)
fit1 <- readRDS("D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Testes\\Usados no artigo\\fit6for.rds")
#fit1 <- bvar.sv.tvp(var1, p=p., tau = tau., nburn = nburn., nrep = nrep., nf = nf., thinfac = thinfac.)
#saveRDS(fit1, "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Testes\\fit6for.rds")
#rm(fit1)

# Primeiro eu monto um objeto só com a previsão
# Eu sei que essa coisa tem três dimensões e quero calcular a média da terceira dimensão e ficar com uma matriz
var6_prev <- apply(fit1$fc.ydraws, c(1,2), mean)

# Agora preciso calcular os erros de previsão em cada período
# nf é o número de períodos de previsão
# amostra1 é onde ficam guardados os valores que estou prevendo
# var1_prev é onde ficam guardados os valores que eu previ

amostra6 <- t(amostra_1)
erro_prev <- amostra6 - var6_prev

q_t <- abs(erro_prev/apply(var1[((tau.+p.+1):nrow(var1)),], 2, FUN = den_MASE))

MASE_6 <- apply(q_t, 1, sum)

MAPE_6 <- vector()
MSE_6  <- vector()

for (i in 1:ncol(var1)){
    MAPE_6[i] <- MAPE(x = var6_prev[i,1:nf.], ref = amostra6[i,1:nf.])
    MSE_6[i] <- MSE(x = var6_prev[i,1:nf.], ref = amostra6[i,1:nf.])
}

Agora estimo um modelo VAR(1) simples para comparar

fit_simples <- VAR(var1, p=p.)
var_simples_prev <- predict(fit_simples, n.ahead = nf.)
vars_prev <- matrix(NA, ncol = ncol(var1), nrow = nf.)

for (i in 1:ncol(var1)){
  vars_prev[,i] <-  (var_simples_prev[[1]][[i]][,1])
}

erro_simples <- amostra_1 - vars_prev

q_t <- abs(erro_simples/apply(var1[((tau.+p.+1):nrow(var1)),], 2, FUN = den_MASE))

MASE_6_simples <- apply(q_t, 2, sum)

MAPE_6_simples <- vector()
MSE_6_simples  <- vector()

for (i in 1:ncol(var1)){
   MAPE_6_simples[i] <- MAPE(x = vars_prev[1:nf.,i], ref = amostra_1[1:nf.,i])
   MSE_6_simples[i] <- MSE(x = vars_prev[1:nf.,i], ref = amostra_1[1:nf.,i])
}

12 passos à frente (1 ano)

rm(fit1)
#Inflação - Produto - Juros - Câmbio - Desigualdade.
var1 <- cbind(capital_trabalho_final*100, ipca, PIB*100, swap, cambio) # The bvar function does not allows data.frames
#var2 <- cbind(capital_trabalho_final, swap, PIB2, cambio, ipca) # The bvar function does not allows data.frames
names(var1) <- c("capital_trabalho", "ipca", "pib", "swap", "cambio")
#names(var1) <- c("ipca", "pib", "swap", "cambio","capital_trabalho") # Ficou uma merda

#var1 <- c(1,2,3,4,5)

pos_ct     <- grep("capital_trabalho", names(var1))
pos_swap   <- grep("swap", names(var1))
pos_pib    <- grep("pib", names(var1))
pos_cambio <- grep("cambio", names(var1))
pos_ipca   <- grep("ipca", names(var1))

nburn.    <- 1000000 # Burn the first 100kk replications
nrep.     <- 100000 # Estimate another 100k times after burn in
thinfac.  <- 4 # But keep only 1 in 4 estimates
p.        <- 1 # Number of lags
tau.      <- 48 - p. # Number of observations to use as prior
nf.       <- 12 # Number of forecast periods
k_B.      <- 6 # Fator de inflação na priori
k_A.      <- 4 # Fator de inflação na priori

# No idea what's this for
tamanho <- floor(1/thinfac.*nrep.)
# número de tempos
tempo <- nrow(var1)

# Prepara um vetor com as últimas nf. observações para deixar de fora
amostra_1 <- var1[(nrow(var1)-nf.+1):(nrow(var1)),] 
var1      <- var1[1:(nrow(var1)-nf.),]

set.seed(2747)
fit1 <- readRDS("D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Testes\\Usados no artigo\\fit12for.rds")
#fit1 <- bvar.sv.tvp(var1, p=p., tau = tau., nburn = nburn., nrep = nrep., nf = nf., thinfac = thinfac.)
#saveRDS(fit1, "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Testes\\fit12for.rds")
#rm(fit1)

# Primeiro eu monto um objeto só com a previsão
# Eu sei que essa coisa tem três dimensões e quero calcular a média da terceira dimensão e ficar com uma matriz
var12_prev <- apply(fit1$fc.ydraws, c(1,2), mean)

# Agora preciso calcular os erros de previsão em cada período
# nf é o número de períodos de previsão
# amostra1 é onde ficam guardados os valores que estou prevendo
# var1_prev é onde ficam guardados os valores que eu previ

amostra12 <- t(amostra_1)
erro_prev <- amostra12 - var12_prev

q_t <- abs(erro_prev/apply(var1[((tau.+p.+1):nrow(var1)),], 2, FUN = den_MASE))

MASE_12 <- apply(q_t, 1, sum)

MAPE_12 <- vector()
MSE_12  <- vector()

for (i in 1:ncol(var1)){
    MAPE_12[i] <- MAPE(x = var12_prev[i,1:nf.], ref = amostra12[i,1:nf.])
    MSE_12[i] <- MSE(x = var12_prev[i,1:nf.], ref = amostra12[i,1:nf.])
}

Agora estimo um modelo VAR(1) simples para comparar

fit_simples <- VAR(var1, p=p.)
var_simples_prev <- predict(fit_simples, n.ahead = nf.)
vars_prev <- matrix(NA, ncol = ncol(var1), nrow = nf.)

for (i in 1:ncol(var1)){
  vars_prev[,i] <-  (var_simples_prev[[1]][[i]][,1])
}

erro_simples <- amostra_1 - vars_prev

q_t <- abs(erro_simples/apply(var1[((tau.+p.+1):nrow(var1)),], 2, FUN = den_MASE))

MASE_12_simples <- apply(q_t, 2, sum)

MAPE_12_simples <- vector()
MSE_12_simples  <- vector()

for (i in 1:ncol(var1)){
   MAPE_12_simples[i] <- MAPE(x = vars_prev[1:nf.,i], ref = amostra_1[1:nf.,i])
   MSE_12_simples[i] <- MSE(x = vars_prev[1:nf.,i], ref = amostra_1[1:nf.,i])
}

24 passos à frente (2 anos)

rm(fit1)
# Previsão
#Inflação - Produto - Juros - Câmbio - Desigualdade.
var1 <- cbind(capital_trabalho_final*100, ipca, PIB*100, swap, cambio) # The bvar function does not allows data.frames
#var2 <- cbind(capital_trabalho_final, swap, PIB2, cambio, ipca) # The bvar function does not allows data.frames
names(var1) <- c("capital_trabalho", "ipca", "pib", "swap", "cambio")
#names(var1) <- c("ipca", "pib", "swap", "cambio","capital_trabalho") # Ficou uma merda

#var1 <- c(1,2,3,4,5)

pos_ct     <- grep("capital_trabalho", names(var1))
pos_swap   <- grep("swap", names(var1))
pos_pib    <- grep("pib", names(var1))
pos_cambio <- grep("cambio", names(var1))
pos_ipca   <- grep("ipca", names(var1))

nburn.    <- 1000000 # Burn the first 100kk replications
nrep.     <- 100000 # Estimate another 100k times after burn in
thinfac.  <- 4 # But keep only 1 in 4 estimates
p.        <- 1 # Number of lags
tau.      <- 48 - p. # Number of observations to use as prior
nf.       <- 24 # Number of forecast periods
k_B.      <- 6 # Fator de inflação na priori
k_A.      <- 4 # Fator de inflação na priori

# No idea what's this for
tamanho <- floor(1/thinfac.*nrep.)
# número de tempos
tempo <- nrow(var1)

# Prepara um vetor com as últimas nf. observações para deixar de fora
amostra_1 <- var1[(nrow(var1)-nf.+1):(nrow(var1)),] 
var1      <- var1[1:(nrow(var1)-nf.),]

set.seed(2747)
fit1 <- readRDS("D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Testes\\Usados no artigo\\fit24for.rds")
#fit1 <- bvar.sv.tvp(var1, p=p., tau = tau., nburn = nburn., nrep = nrep., nf = nf., thinfac = thinfac.)
#saveRDS(fit1, "D:\\Onedrive - Aisha\\OneDrive\\Documentos\\Mestrado Economia\\Dissertação\\YSI - Buenos Aires\\Testes\\fit24for.rds")
#rm(fit1)

# Primeiro eu monto um objeto só com a previsão
# Eu sei que essa coisa tem três dimensões e quero calcular a média da terceira dimensão e ficar com uma matriz
var24_prev <- apply(fit1$fc.ydraws, c(1,2), mean)

# Agora preciso calcular os erros de previsão em cada período
# nf é o número de períodos de previsão
# amostra1 é onde ficam guardados os valores que estou prevendo
# var1_prev é onde ficam guardados os valores que eu previ

amostra24 <- t(amostra_1)
erro_prev <- amostra24 - var24_prev

q_t <- abs(erro_prev/apply(var1[((tau.+p.+1):nrow(var1)),], 2, FUN = den_MASE))

MASE_24 <- apply(q_t, 1, sum)

MAPE_24 <- vector()
MSE_24  <- vector()

for (i in 1:ncol(var1)){
    MAPE_24[i] <- MAPE(x = var24_prev[i,1:nf.], ref = amostra24[i,1:nf.])
    MSE_24[i] <- MSE(x = var24_prev[i,1:nf.], ref = amostra24[i,1:nf.])
}

Agora estimo um modelo VAR(1) simples para comparar

fit_simples <- VAR(var1, p=p.)
var_simples_prev <- predict(fit_simples, n.ahead = nf.)
vars_prev <- matrix(NA, ncol = ncol(var1), nrow = nf.)

for (i in 1:ncol(var1)){
  vars_prev[,i] <-  (var_simples_prev[[1]][[i]][,1])
}

erro_simples <- amostra_1 - vars_prev

q_t <- abs(erro_simples/apply(var1[((tau.+p.+1):nrow(var1)),], 2, FUN = den_MASE))

MASE_24_simples <- apply(q_t, 2, sum)

MAPE_24_simples <- vector()
MSE_24_simples  <- vector()

for (i in 1:ncol(var1)){
   MAPE_24_simples[i] <- MAPE(x = vars_prev[1:nf.,i], ref = amostra_1[1:nf.,i])
   MSE_24_simples[i] <- MSE(x = vars_prev[1:nf.,i], ref = amostra_1[1:nf.,i])
}

Agora vamos juntar todo mundo numa tabelinha bonitinha.

	capital.trabalho	ipca	pib	swap	cambio
MASE TVP-VAR 1	0.0049	0.4644	0.4324	0.1455	0.6737
MASE VAR 1	0.3216	0.6578	0.3802	1.7425	0.5887
MSE TVP-VAR 1	0.0020	0.0261	1.0462	0.0048	3.7993
MSE VAR 1	8.8195	0.0523	0.8089	0.6843	2.9015
MAPE TVP-VAR 1	0.0009	0.0568	4.3178	0.0106	0.0191
MAPE VAR 1	0.0596	0.0804	3.7968	0.1265	0.0167
MASE TVP-VAR 6	1.6333	4.1112	3.1941	2.6224	6.5894
MASE VAR 6	10.7608	3.3711	7.1403	24.1303	4.7334
MSE TVP-VAR 6	7.5701	0.0918	1.8476	0.0748	14.0991
MSE VAR 6	3.1085	0.1935	1.8125	10.2476	1.2028
MAPE TVP-VAR 6	0.0502	0.0857	1.5908	0.0313	0.0320
MAPE VAR 6	0.0300	0.1329	1.5482	0.4413	0.0091
MASE TVP-VAR 12	5.2766	63.7988	5.9903	83.3788	30.8935
MASE VAR 12	33.2865	27.9069	16.8345	71.1768	56.7461
MSE TVP-VAR 12	23.8266	3.9851	1.9681	12.9829	69.9769
MSE VAR 12	8.8536	4.1032	1.9578	31.1229	17.9026
MAPE TVP-VAR 12	0.0801	0.6638	1.5928	0.4450	0.0774
MAPE VAR 12	0.0498	0.6726	1.5042	0.6905	0.0387
MASE TVP-VAR 24	141.4759	308.2861	32.3579	246.3455	107.3844
MASE VAR 24	96.6396	115.2667	40.2186	102.1564	410.8116
MSE TVP-VAR 24	10957.0877	25.3176	28.0481	54.9363	199.3617
MSE VAR 24	39.2327	22.8249	2.4188	26.1457	227.8310
MAPE TVP-VAR 24	1.1113	1.2500	6.1358	0.6622	0.1353
MAPE VAR 24	0.0642	1.2224	1.1793	0.4946	0.1474

Artigo YSIS - Política Monetária x Desigualdade

Aishameriane Schmidt

Maio 2018

Dúvidas pra ver com o Guilherme

Introduction

Downloading packages

Data handling

Downloading data

Downloading the data

Visualizing data

Treating for seasonality and structural breaks

1. Capital-Labor Ratio

2. Selic (interest rate)

2. Swap DI pré 3 meses

3. PIB mensal

4. Exchange rate

5. IPCA

Variávels finais após ajuste sazonal

Descritivas

Estimating the BVAR

Funções impulso resposta

Impactos dos choques monetários

Em todas variáveis do modelo, no último período

Na razão capital trabalho

No PIB

No Câmbio

No IPCA

Choques do PIB

Na razão capital trabalho

No Câmbio

No swap

No IPCA

Choques na taxa de câmbio

Efeito do câmbio na razão capital-trabalho

Efeito do câmbio No Swap

Efeito do câmbio No IPCA

Efeito do câmbio No PIB

Choques no IPCA

Efeito do Ipca na razão capital-trabalho

Efeito do Ipca no swap

Efeito do IPCA na taxa de câmbio

Efeito do IPCA no PIB

Avaliando os parâmetros

Razão Capital-Trabalho

Swap

PIB

Taxa de Câmbio

IPCA

Plotando as variáveis todas juntas

Gráfico dos coeficientes das outras variáveis defasadas na equação da razão capital-trabalho

Diagnóstico

Diagnósticos para as médias

Diagnóstico da autocorrelação

Diagnóstico multivariado

Previsão

Medidas de erro de previsão usadas

1 passo à frente

6 passos à frente (1 semestre)

12 passos à frente (1 ano)

24 passos à frente (2 anos)